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Abstract. We present a formalism to describe slowly decaying systems in the context 
of finite Markov chains obeying detailed balance. We show that phase space can 
be partitioned into approximately decoupled regions, in which one may introduce 
restricted Markov chains which are close to the original process but do not leave these 
regions. Within this context, we identify the conditions under which the decaying 
system can be considered to be in a metastable state. Furthermore, we show that 
such metastable states can be described in thermodynamic terms and define their free 
energy. This is accomplished showing that the probability distribution describing the 
metastable state is indeed proportional to the equilibrium distribution, as is commonly 
assumed. We test the formalism numerically in the case of the two-dimensional kinetic 
Ising model, using the Wang-Landau algorithm to show this proportionality explicitly, 
and confirm that the proportionality constant is as derived in the theory. Finally, 
we extend the formalism to situations in which a system can have several metastable 
states. 
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1. Introduction 

The description of macroscopic states has been achieved succesfully in the case of systems 
which are in thermodynamic equihbrium. Indeed, for these, Gibbs' approach via the 
canonical or other ensembles, describes in a well-defined manner most equilibrium 
systems (under appropriate ergodicity assumptions), and is also well-supported by 
experiment. On the other hand, the goal of similarly describing appropriate macroscopic 
states for systems out of equilibrium has, so far, not been achieved in general, and the 
attempts to do so have met with serious difficulties. 

An interesting intermediate case is found for so-called metastable systems. These 
arise typically near first order phase transitions, when a phase, which is not the 
thermodynamically most stable one for the values of the parameters considered, is 
nevertheless observed to persist over very long times. A typical instance occurs, for 
example, when a system is rapidly quenched from a region in which a given phase is 
stable to another region in which the phase is unstable. 

Traditionally (see, for example. Maxwell PP), such states have been considered 
in a purely thermodynamic framework. Analytic forms of the free energy, which 
involved non-convexities, were taken to refiect physical reality, and the states showing 
local thermodynamic stability (for example, those which have positive susceptibility 
or compressibility) were identified with metastable states as observed in nature. This 
straightforward explanation, however, suffered a serious blow when it was understood 
that short-range systems must necessarily have convex free energies, as was shown, for 
example, in The presence of non-convex parts in the free energy of a van der Waals 
system is hence necessarily linked to the fact that the van der Waals model is exact only 
in the limit of long-range potentials. 

Nevertheless, a great deal of work followed along these lines. In particular, Langer 
|3j showed how the original program of Maxwell and van der Waals could be carried 
out under quite reasonable assumptions on the nature of the singularities present near 
the coexistence curve of two phases. The metastable phase was there assumed to be 
destabilized by the possible - but very unlikely - presence of large droplets of the stable 
phase. That droplets of the stable phase are unlikely to arise was explained by the 
fact that they must appear previously as smaller droplets, which are themselves highly 
unlikely, since they are strongly suppressed by the effect of surface tension. In such 
models, the free energy and the partition function are obtained from the corresponding 
equilibrium quantities by a highly non-trivial analytic continuation around an essential 
singularity. The reason why the physical metastable free energy should be identified 
with the result of such an analytical continuation, however, is not entirely clear in this 
approach. Indeed, this could hardly be otherwise, since the basic objects with which 
this framework operates are drawn from equilibrium statistical mechanics, whereas 
the metastable state, being subject to decay to a state very different from itself, is 
intrinsically a non-equilibrium object, albeit perhaps a particularly simple one. 

In a different line of attack, Penrose and Lebowitz |1] studied metastability directly 
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from a dynamical point of view. Their considerations were limited to rather particular 
models, but the principles involved seem capable of great generalization. The crucial 
idea was the introduction of a restricted ensemble, defined in such a way as to prevent 
nucleation from ever taking place. If such a dynamics could indeed be defined in 
such a way as to have the characteristics of a typical equilibrium process, then the 
thermodynamic approach would be, in a sense, vindicated. 

In more recent work, Gaveau and Schulman [S] have succeeded in making this quite 
precise in the very general framework of arbitrary Markov processes. Their approach 
consists in assuming that some non-equilibrium eigenstate of the linear operator arising 
in the master equation describing the dynamics has an anomalously small eigenvalue 
with respect to all others. In systems for which the thermodynamic limit has not yet 
been taken, arguments taken from the classical theory of nucleation will often strongly 
suggest the existence of such slowly decaying states. As to the thermodynamic hmit, it 
has been forcefully argued in [3] that it may in fact be quite misleading when applied 
to metastable systems. In their work, they essentially show the following: 

(i) Under their assumptions, the space of configurations separates into two disjoint 
subsets, which are both almost invariant under the dynamics, one of which can be 
identified as the metastable region. In particular, if the system is started in this 
region and allowed to evolve until it attempts to leave it, after which it is killed, 
then the evolution of the system is similar to that in which the system evolves in 
the usual manner. 

(ii) Similarly, they show that the average time to escape from the metastable region is 
very large. 

In this work we develop a rather simple formalism, along lines similar to those of 
Gaveau and Schulman, that sharpens and extends their results in various ways; our 
development was outlined in 0. We confine ourselves to ergodic and acyclic processes 
satisfying detailed balance, for which the stationary distribution can unambiguously be 
identified with the equilibirum distribution in the appropriate ensemble of statistical 
mechanics. For such systems, a restricted dynamics in which the process is refiected, 
instead of killed, each time it attempts to leave the metastable region is shown to also 
be close to the original process. (This kind of restriction was also used in a slightly more 
specific context in This restricted process reaches an equilibrium state described 
by a distribution which is proportional to the equilibrium distribution of the original 
unrestricted process. We also show that under suitable conditions, any initial condition 
decays quickly into either a metastable state, in which case the system is described by 
the equilibrium distribution of the restricted process, or to equilibrium. 

In this respect, another important point in which we sharpen the results of Gaveau 
and Schulman is the following: they prove that on average it takes a long time for a 
system to decay from the metastable phase to equilibrium. This, however, is not enough 
to account for the expected behavior of metastable systems: for example, if a system 
were to decay in a time of order one with probability 1 — e and in a time of order 
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with probability e, then the average decay time would be of order e~^, and hence 
would diverge as e — > 0, yet no one would call such a system metastable. Indeed, in 
agreement with [Zj, one expects metastable states to reach a (quasi-)stationary regime 
quickly and then, in a relatively abrupt maner, decay to equilibrium. The broadness of 
the distribution of the time intervals during which the state remains in the metastable 
(quasi-) stationary state is what gives rise to the slow decay of the distribution describing 
such a state. Here we show that, under the assumptions required to define metastability, 
it is in fact true that the probability of decay in a short time is very small, in accordance 
with the intuitive picture given above. 

These result support the idea behind the restricted process approach. Furthermore, 
since we are really dealing with partition functions of two systems, both of which are 
equilibrium systems, the logarithm of the proportionality constant is related to the 
difference in free energy between the unrestricted and restricted systems, corresponding 
to the stable and the metastable phases respectively. 

Another interesting consequence of these observations is the following: since a 
metastable system can, to a good approximation, be described by an equilibrium 
process over certain time scales and the usual connections between time correlations 
and response to small external perturbations (fluctuation-dissipation theorem) hold 
exactly in the restricted dynamics, then, again to a good approximation, they will also 
hold in the metastable state. It is therefore legitimate, say, to measure the frequency 
dependent susceptibility in a metastable state by computing the Fourier transform of 
the magnetization autocorrelation function 

To illustrate some of the results described above, we study a two-dimensional Ising 
model subject to an external field. We parametrize the phase space by reduced variables 
(in this case magnetization and total spin-spin interaction energy, which are adequate for 
the system sizes we are considering) and evaluate the equilibrium distribution over the 
complete parameter space using the Wang-Landau algorithm. Within the metastable 
region, we compare the equilibrium distribution to the metastable distribution obtained 
from Monte Carlo simulations of the kinetic Ising system. In the metastable region 
within the space of reduced variables, we show that the metastable distribution is indeed 
proportional to the equilibrium distribution, with the proportionality constant being as 
derived in the theory. 

Finally, we extend the formalism to the case in which the system has several 
metastable states. This gives rise to minor complications due to the possibility that 
the system may decay to equilibrium by passing through other metastable states. 

The outline of the paper is as follows. In Section |2l we review the formalism, which 
is similar to that used by Gaveau and Schulman, and the assumptions and notation to 
be used throughout the paper. In Sectional we present and derive the results described 
above in the case where the system has a single metastable state. In Section |31 we show 
how our ideas can be applied to the Ising model, at least for sufficiently small systems. 
The results shown in the previous sections can be confirmed using this test model. In 
Section we discuss the complications appearing when a finite number of metastable 
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states are taken into account, instead of only one. Finally, in Section IHl we present some 
conclusions and outlook. 

2. Theoretical framework 

We set up the description of slowly decaying as well as metastable states within the 
general framework of Markov process, which can then be applied to a large variety of 
systems. Given a finite set F of elements a, we consider a continuous time Markov chain 
on this set defined by transition rates Wa->a'- The probability P{<7,t) to encounter the 
system at time t in the configuration a then obeys the master equation 
dP 

— {a, t) = Yl [^-'--^(^'' t) - W^^^^P{a, t)] . (1) 

If this Markov process satisfies the conditions of ergodicity and aperiodicity, see 0, 
which are usually satisfied in the systems we are interested in |, then the probability 
distribution P{a,t) approaches a unique equilibrium distribution PQ^a) as t — oo. 
We will further assume that the system obeys detailed balance. That is, 

W„,^M(^') = W„^,.PM (2) 
holds for all a. We rewrite ^ in the operator form 
dP 

where P is a vector with index a, and L is a linear operator on the space of all such 
vectors. A scalar product of two vectors and ip can be defined as 

(0,^):=L-^' (4) 

under which, given the detailed balance condition (|21), the operator L is self-adjoint: 

(0,L^) = (L0,^). (5) 

Since the underlying vector space is finite-dimensional, there is a complete 
orthonormal set of eigenvectors Pn satisfying 

LPn = -n^Pn, (6) 

where the fin are by definition arranged in increasing order, and N is the number of 
elements of F, i.e. the number of possible configurations. The existence of an equilibrium 
distribution implies that fio = and the corresponding Pq is in fact the equilibrium 
distribution. All other f2„ are strictly positive. 

Using the orthonormality of the P„ we can write 

(Po,P„) = 5^P„(a)=5„,o, (7) 

cr 

X Note that we are dealing with finite systems only, so that problems of ergodicity due to, say, phase 
transitions do not arise. 
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implying that Po{(j) is normalized and that adding to it arbitrary multiples of Pni<^), 
when n > 1, does not alter this normalization. The completeness of the eigenvectors 
(jni) implies that 

(5^„((t) ■=6^,,^= y — — . (8) 

This leads to an exact expression for the probability of arriving from do to cr in time t: 
P(a, t; ao, 0) = e^X, (a) = Po(a) + ^iMf^e'^^t^ 

In the following, we shall say that a system is slowly decaying if at least one of its 
eigenvalues f2„ is much less than all the others. At first, we shall limit ourselves to the 
case in which there is only one slow eigenvalue, that is, when Qi <^ Qn for all n > 2. 

Now consider a process evolving from the initial condition a^. Then, from ©, in 
the relevant time range il^^ -C t -C ^i^, one finds that the configuration a is occupied 
with the following (time-independent) probability 

P{a)=Po{a) + ^^P,{a). (10) 

Note that, due to ((Zj), this is normalized. Also, since it differs exponentially little from 
the exact result, we may conclude that it is positive, except perhaps in some places 
where it assumes exponentially small negative values. This situation can be corrected 
by setting the negative values to zero and recomputing the normalization, which leads 
to negligible alterations. 

This result focuses our attention on the value Pi(cro)/Po(cro), which characterizes 
the nature of the initial condition. This quantity will be central to all that follows. 
In particular it will allow us to determine when the initial condition can be called 
metastable and the resulting probability distribution given by (fTUj) can justifiably be 
identified with that of a metastable state. 

In what follows, we denote Pi(cr)/Po(cr) by C(cr), and the maximum value of C(cr) 
over all 0" G r by C. Next we define the sets Fm and Feq as follows: 

T^:=<aeT:^<^<C'>, (111 



2 - Po(^) 

and Feq is defined as the complement of Fm in F. We show in Section IHl that the choice 
of the factor 1/2 to define the lower bound on C{a) in pi|) is relatively arbitrary. 

Equation (fTTjl defines a partition of phase space into two disjoint sets. In 
the following we shall address the question of the extent to which we can use this 
partition to define a metastable state, and in particular, to understand when a standard 
thermodynamic approach to the study of such systems is legitimate. 

To this end we single out among the slowly decaying systems those characterized by 
the condition that the probability of being found within Fm in equilibrium is negligibly 
small, i.e. such that 
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We will call metastable systems the slowly decaying systems satisfying this condition. 
In particular, from normalization, metastable states satisfy: 

5^ PoM = 1-/^^1. (13) 

We now turn to proving various properties both for slowly decaying systems not 
satisfying condition (fT^ . and for metastable systems, with a view to justifying the 
usual assumptions concerning the description of the latter. 



3. Results and proofs 

We begin by considering a slowly decaying system with a single slow mode, so that its 
phase space is partitioned into Feq and Fm as before. We first consider the case in which 
the initial condition cxo satisfies 

C(ao) = C. (14) 

In this case, we define Qi (a) as the quasi-stationary distribution which arises from this 
initial condition over the time range fig ^ -C t ^ ^i^, given by (see (jlUj) ) 

Q,{a) ■.= Po{a) + CP,{a). (15) 



3.1. Probability of exit from the metastable state 

Let us define the random variable T as the time at which a path starting at ao satisfying 
()14|1 reaches Feq for the first time. A key result is that with high probability this time 
is large. Indeed, for these processes, we have 

P(r < t) < 2 (l - e"^^*) = O(fiit). (16) 

Here P(- ■ ■) denotes the probability of an event; for example, the LHS of ()16p denotes 
the probability of T being less than t. 

To prove this, we proceed as follows. We denote by a(t) the path followed by the 
process starting at a{0) = ctq of the Markov process defined by (H)), and by E{- • ■} the 
expectation value. 

The set is defined by a condition on the function C(cr), so to study the first exit 
time T from this set, we must consider the evolution of C[cr(t)] as a function of time. 
We thus consider, for t' > t, 

E {C[a{t')] I a{t)} = J2 C{a')P{a', t'\a, t) 



= e^i(*-*')C[(T(t)], (17) 

where we have used equation (jU)) and the definition of C{a). The last equality follows 
from the orthonormality of the basis P„(cr). Thus we have 

E|e^^*'C[a(t')]|ff(t)} = e^^*C[a(t)], (18) 
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so that e^^*C[cr(t)] is a martingale jU]; intuitively, this means that, on average, it neither 
grows nor decreases with time. 

Furthermore, T is a so-called stopping time, that is, it is known at time t whether 
the event T < t has taken place or not. If we now define r = min(t, T), then by standard 
theorems on martingales and stopping times (see e.g. jH]), it follows that 

E{e^^-C[a(r)]} =C(ao), (19) 

where (Tq is the initial condition of the process, which we chose such that C(cro) = C, 
its maximum possible value. We can therefore estimate the LHS of (|T9|) from above: 

C = E {e^^"C[(T(r)] } < { e^'^\ T < t} F{T < t) +E { e^''C[a{t)] \T>t}F{T> t) 



< Ce'''' I -P(T <t) + [l- P(T < t)] I , (20) 

from which ()16|) follows immediately. 

This result is of considerable interest. It represents a significant sharpening of a 
result due to Gaveau and Schulman stating that the average value (T) is large. Here 
we show that, for appropriate initial conditions, the system is very unlikely to leave Fm 
before time t in the relevant time range t <C From this result one may also derive 
the following estimate on Poic) and Piic), which will be of use later: 

iy:= Yl ^i(^) = Yl l^o(^) + ^ < t) < 1. (21) 

The inequality follows from the fact that u is equal to the total probability of finding a 
system started at an initial condition ctq with C(o"o) = C in Fgq at time t. But this is 
less than P(T <t), so that the estimate (PT|) follows. If we think of Qi(cr) as describing 
a metastable state, then ^IT\ states the (perhaps unsurprising) fact that the metastable 
state is entirely concentrated outside Feq. Note that the converse, namely that Pq has 
only negligible weight in Fm cannot be shown in a similar way. Rather, this condition 
is what we introduce in equation (fT^ as an additional hypothesis to single out true 
metastable states from slowly decaying systems. 

3.2. Definition of restricted process in the metastable state 

We now introduce a restricted Markov process in order to be able to treat the slowly 
decaying system as if it were in fact in equilibrium. To this end, define the following 
restricted transition rates: 

yyR j Wa'^a cr, a' e Fin or (J, a' E Feq ^22) 

"'^^ ■ I otherwise. ^ ' 

It is clear that the rates W^_^^ only allow for connections within F^ or Feq. In fact, the 
R process can be intuitively understood as a process that imposes refiecting boundary 
conditions at the border separating from Feq §. Since Po{cr) satisfies detailed balance 

§ This process differs from the one considered in [H] , in which the process is killed whenever it attempts 
to leave Fn,. 
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in the original process, it is still the equilibrium distribution for this restricted process, 
with respect to which it also satisfies detailed balance. But the restricted system is no 
longer ergodic, and therefore Po{<j) is not the unique stationary distribution. Indeed, 
Pf{a) defined by 

y b io(crj, G 1 cq 

is stationary for any constants C and S' . In particular, we may choose these constants 
so that Pi{(^) = and (P/^, P^) = 1. This implies that 

c = ^^^vfr ; = - V^'- (24) 

\22r^Po{(^) J 

Of course, it is now very tempting to identify P/^ with Pi. In order to do this, we 
need to show that the process defined by (j22j) . which we denote by R (for Restricted), 
remains close to the original Markov process defined by the rate W^^„' , which we denote 
by P (for Physical). This can indeed be shown for a slowly decaying system, if the initial 
condition ctq satisfies C{ao) = C and t is in the relevant time range Qit <C 1 <C Q2t- We 
define closeness as follows: for any subset X C Fm, define 

px{t) := |P{ap(t) G X} - P W(t) G X}| , (25) 

where ap{t) and aR{t) are paths of the P and R processes, respectively. We will say 
that the two processes are close in variation if px{t) is small for any X C Fm. 

For the proof, we make the following observation, inspired by the coupling 
techniques of probability theory. We define a compound process K = (ap, ap) on 
the product space F x F as follows: ap moves according to the P process, that is, via 
the rates W, and ap follows ap around as long as the latter remains in Fm. As soon as 
ap leaves T^, however, each process evolves independently according to their respective 
rates. By construction, the projections of the compound process K on either subspace 
yield the processes R and P, respectively. The two paths ap{t) and ap(t) can thus be 
viewed as projections of the process K. 

We wish to show that sup px{t) is small in the relevant time range. This is achieved 

xcr 

as follows: again let us introduce the random time T as the first time at which ap(t), 
starting from uo for which C{ao) = C, leaves Fm. It then follows by the construction of 
the process K that 

p^{t) = |P {ap{t) eX\T<t} + ¥ {ap{t) eX\T>t} 

- F{ap{t)eX\T<t}-F{aR{t)eX\T>t}\ 

= \f{ap{t) e X\T <t} ~¥{aR{t) eX\T <t}\ 

= \F{ap{t) G X\T < t}-F{aR{t) G X\T < t}\ x P(T < t) 

< P(r < t), (26) 

which is indeed small for t <^ according to (jlfij) . As this holds for any X C F, 
the probability distribution for the restricted process is close in variation to that of the 
physical process. 
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This implies that within the relevant time range, if ctq G Fm, then 

Pp(a,t;ao,0) ^Pfi(a,t;ao,0), (27) 

where Pp,i?((T, t; do, 0) denote the transition probabilities from do to a in a time t 
for the physical and the restricted process respectively. Again choosing do such that 
C(o"o) = C and expressing each of these distributions in terms of the eigenfunctions of 
their respective evolution operators, we have: 

Pp(a, t- ao, 0) = + CPf (a)e-^f * + 5£MSMe-^.-* (28) 



and 



n=2 



Pp(a, t; ao, 0) = P,{a) + C'P/^(a)e-^?* + El^^}El^e-^'-K (29) 

Then, the relation expressed in (j?rj) implies that P^io-) ~ Pni^)^ these quantities 
are not negligible in T^, in which case we also have ~ Q^. Thus, again in the time 
range fig ^ <^ ^i^, we are left with 

gf (a) := Po(a) + CP[{a) ^ Po(a) + C"P«(a) =: Qf (a), (30) 
which, together with the fact that Xlrm Qii^) ~ ('-'^) ~ l^ads to 

C^C and Pf{a) ^ P^{a). (31) 

We have therefore two results of interest: on the one hand, the first passage time from 
a state ctq satisfying C(cro) = C to Feq is very unlikely to be short. On the other, the 
process starting at ctq restricted to remain forever in is quite similar to the original 
unrestricted process for times in the relevant time range. 



3.3. Generalisation to other initial conditions 

So far we have restricted attention to initial conditions ctq such that C(o"o) = C. We now 
show that to a large extent this requirement on the initial condition becomes unnecessary 
if one makes the hypothesis that the system is a metastable one, that is, that (jl2j) 
holds. For such systems we will prove the following basic property: no matter what 
the initial condition ctq is, provided it satisfies C{crQ)/C = 0(1), within a time of order 
0,2^ the system will either be in Feq or else it will satisfy approximately the condition 
C[cr(t)] = C. This means, therefore, that the two results described above can be applied 
whatever the initial condition, provided only that the process remains within for a 
short time. 

To show this we first need an auxiliary result, which also depends on the extra 
assumption that defines metastability: If we consider an initial condition ctq^^ such that 
C(o"g^^) = (1 —p)C, then the probability that this initial condition ends up in Feq, after 
a time significantly larger than fig ^ has elapsed, is p. Indeed, in the relevant time range 
f2^^ ^ t ^ ^2 ^5 this probability is 

F{a^^\t) e Feq|C[4^)] = {1-P)C] ^ (1 -p) ^ Qiicj)+p J2 Po{<y)^P, (32) 
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where we have combined the facts that for slowly decaying systems Qi has essentially 
no weight in Fgq (equation (j7T|) ). and that Pq has no weight in Fm (equation (HH)). 

This result further implies that for values of p such that p/u ^ 1, where u is defined 
by (EU, we have 

F(p):= Yl (33) 

<t;C((t)<(1-p)C 

Indeed, consider a system evolving from an initial state given by P{(J, 0) = Qi{(j). The 
probability of finding the system in Feq after a time t will then be given by 

J2 P(a,t) = ^P(a(t)GFeq|ao = a)P(a,0). (34) 

(Tgrcq crgr 

It then follows that 

J2 ^(^' ^) ^ ^ W ^ Tcqko = ^) i^(^, 0). (35) 

o-eFeq <t;C((t)<(1-p)C 

Now, equation (jH^ implies that, in the relevant time range, 

P(a(t)GFeq|ao = a)>p, if C{a) < {1 - p)C, (36) 

so that 

5^P(a,t)>p Yl ^(^'0)- (37) 

O-ercq <T:C(cr)<(l-p)C 

However, since the initial state was Qi, which is essentially stationary in this time range, 
we have P{<7,t) ^ Qi{<y), giving 

pFip)< J2 Qii^) = ^- (38) 
o-ercq 

Since u is negligibly small, we thus find that F{p) < v/p ^ 0. 

In a similar way, we can show that -Po(o') is non- negligible only for the states for 
which C{a) ^ 0. This time consider 

Gip):= ^o(^)- (39) 

After a time of order \ ^it least (1 —p)G{j>) of these states will end up in Fm- Thus, 
following the same line of reasoning as before, we can conclude that 

{l-p)G{p)< 5^Po(^)=/x. (40) 

But our basic hypothesis is that for metastable states, is negligibly small, thus 
G{p) ^ 1 for p ^ 1 — /i. In other words, if the condition for metastability, equation (fT^ . 
holds when F^ is defined by the inequalities then a similar claim can be shown 
when the prefactor 1/2 is replaced by essentially any other number well within v and 
1 — yU. This means, in fact, that outside a boundary set with relatively small measure 
both with respect to Pq and to Qi) the function C{a) takes only the values and C . 

Conversely, it is obvious that if G{p) -C 1 for all p within u and 1 — /x, then the 
state will be metastable in the sense that equation (jl3p is satisfied. No similar converse 
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statement holds for F{p): in that case, F{p) was found to be neghgible as a consequence 
of the fact that the probabihty of leaving within the relevant time range is negligibly 
small. This, as we have seen, is the case for arbitrary slowly decaying systems, if the 
initial condition do satisfies C{ao) = C. Thus, for non-metastable slowly decaying 
states, F{p) would be negligible due to the very slowness of their decay. However, such 
modes would not correspond to metastable states unless assumption (fT^ held, and thus, 
G{p) < 1. 

For non-metastable slowly decaying systems, we have that G{p) is not negligible, 
indicating that one can find states a in equilibrium with any value of C{a), including 
C{a) ~ C, all of which would decay slowly. For such initial conditions, the almost 
certain absence of decay within the relevant time range cannot be expected. In fact, 
instead of reaching a stationary state which suddenly decays, the properties of these 
systems will evolve continuously in time until they reach equilibrium. Physical examples 
of such slowly decaying systems are hard to come by: the obvious instances that come 
to mind (slow hydrodynamic modes, such as diffusion, necessary to reach a uniform 
equilibrium from a long- wavelength perturbation) almost invariably involve a quasi- 
continuum gapless spectrum near zero, and are thus ruled out by our basic assumption. 
On the other hand, a trivial, though unenlightening, example shows that non-metastable 
but slowly decaying states do exist: if one specific spin in an Ising model is flipped at 
a much slower rate than all others, it will, as is easily verified, create a slowly decaying 
eigenstate which is not metastable in the sense that it does not satisfy ()12|) 



3.4- Structure of metastable states 

The picture that emerges then, is that for systems having metastable states, after a 
relatively short transient time, the system will only be found in states a for which 
either C{a) ~ (equihbrium) or C{a) ~ C (metastable), independently of the initial 
condition. Further, the dynamical behavior is described to a good approximation by 
the restricted Markov process which is reflected whenever it attempts to go from Fm to 

r 

eq- 

Finally this can be interpreted as follows: the state in which a metastable state 
remains throughout the relevant time range Qit ^ 1 ^ is determined by Qi, which 
is in principle defined in a way that depends on the dynamics. However, as we have 
seen, it turns out that 

P,{a) = CPo{a) for a G F,,, (41) 
from which immediately follows 

Qiia) = Z.Poia) for a G F^, (42) 
where Zi = l+ C^. 

It follows that the only influence of dynamics on the metastable state is that it 
defines the extent of Fm. In other words, the equihbrium ensemble restricted to a 
suitable subset F^ of phase space describes the metastable state. From this follows, in 
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particular, that one may straightforwardly define thermodynamic quantities such as the 
partition function by 

where the last equality follows from the normalization of Qi and equation ()42j) . as well 
as the fact that, as follows from (fT^ . the term X^o-er -^o(c) is in fact negligible. Note 
that this implies in particular that C ^ 1. 

Similarly, we can show that the fluctuation-dissipation theorem, see for example 
[S], will hold for metastable states, since the dynamical correlation functions over the 
relevant time range will be described by a Markov process close to the one that reflects 
the system back to whenever it attempts to leave it. This process, however, is a well- 
deflned Markov process satisfying detailed balance which has the normalized restriction 
of Po{cr) to Fm as an equilibrium state, so that the fluctuation-dissipation theorem can 
be shown for it in a straightforward way. 



4. An illustration: The kinetic Ising model 

We now proceed to show how these ideas can be applied concretely in the case of the 
two-dimensional Ising model. Here, the conflgurations are collections a = (crj)j=i_...,Ar of 
spins (Tj = ±1 at site i, with energy given by the Hamiltonian 

n{a) = - ^ ataj - h^ai =: E{ct) - hM{a), (44) 

where the flrst sum is over nearest neighbours in an := L x L square lattice with 
periodic boundary conditions, and h is the external magnetic fleld. E{(t) and M{a) are, 
respectively, the spin-spin interaction energy and the magnetization of the conflguration 
a. 

To obtain a kinetic model which can exhibit metastability, we must impose a 
dynamics on the system. For concreteness we use discrete-time Metropolis spin-flip 
dynamics (TU]: spin flips are proposed at random, and accepted with probability 
min{l, exp(— /3Ai7)}, where Aif is the change in the Hamiltonian ()44p due to the 
flip, and j3 := 1/T is the inverse temperature fO]- Note, however, that the only thing 
expected to change under a different local || dynamic rule is the extent of the metastable 
region. The Metropolis dynamics gives a discrete-time Markov chain with a unique 
equilibrium distribution at flxed T given by the canonical distribution 

P,[a) = ^eM-m{cr)i (45) 

where 

Z:= J]exp[-/3H(a)] (46) 

a 

II This caveat is necessary since certain non-local dynamics for the Ising model, such as the Swendsen- 
Wang algorithm JO], suppress metastability altogether. 
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is the partition function, from which we can obtain all thermodynamic information at 
equilibrium. The Hamiltonian (03)), together with such a spin- flip dynamics, gives the 
kinetic (or stochastic) Ising model. 

As is well known, if we fix a subcritical temperature T < T^, and a weak external 
magnetic field h is applied, taken negative (downwards) without loss of generality, 
then the spontaneous magnetization in equilibrium points in the direction of that 
field. However, if we initialize the system with all spins up, then for a broad range 
of parameters, the system remains in this thermodynamically unfavorable positively 
magnetised state for a given (random) length of time, whose mean depends on the 
temperature T and the external field h fl] . This state is the prototype of the metastable 
states we aim to describe. 

Since the Metropolis Markov chain is ergodic and acyclic jHj, the formalism 
developed in the previous sections (when rewritten for discrete-time systems) applies to 
this system. Intuitively, it is clear that the kinetic Ising model started in the metastable 
region has a hierarchy of relaxation times, with one (the escape time from the metastable 
region) being much longer than the others. Assuming that this is reflected in the 
spectral properties required in the derivations above, in this section we show that the 
formalism indeed provides a good description of this metastable state. We remark that 
many rigorous results have been proved on metastability in the Ising model in the low- 
temperature limit: see for a comprehensive review; in particular, the separation of 
eigenvalues required in our formalism has been proved in this limit in fn\ IT^ . However, 
our formalism is valid for any temperature, provided that the eigenvalue separation is 
satisfied. 

4.1. Reduced phase space 

To obtain confirmation of the results of Section IHl in the case of the kinetic Ising model, 
we must identify the metastable and equilibrium regions Fm and Feq and compare the 
equilibrium and metastable distributions in each of these regions. However, given the 
huge size of phase space even for small systems, this program cannot be carried out. 
Instead we must resort to a reduced description of the complete phase space in terms of 
a few variables which, if accurate enough, will reflect the relations we predict over the 
complete phase space. 

Due to the numerical techniques we use (discussed below), we are restricted to 
studying relatively small systems. For ferromagnetic Ising models of such sizes, it 
follows from elementary nucleation theory that E and M are sufficient to characterize 
Fm. Indeed, we know that nucleation occurs whenever a droplet of approximate size 
Rc{P, h) arises spontaneously, where Rc depends on (3 and h, but not on the size of the 
system. Inside the critical droplet, the magnetisation has approximately its equilibrium 
value, whereas outside it has the (generally quite different) metastable value. For small 
systems, it is therefore generally impossible for a critical droplet to appear without 
significantly modifying the magnetisation M. For larger systems, it would be necessary 
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to restrict not only M, but also all magnetisations restricted to cells of size of order Re, 
such restrictions presumably define T^. This has been treated in detail in particular in 
|3]. In the following, since we are limited to small systems, we decribe Fm entirely in 
terms of E and M. 

We refer to the set {a G F : E{a) = E; M{a) = M} of configurations with given 
values of the macroscopic variables E and M as the {E, M) macrostate. In this section 
we work exclusively on a coarse-grained level in terms of such macrostates, for the 
reasons just described, by summing over all configurations a belonging to a macrostate. 
For example, summing over the {E, M) macrostate, we obtain 



is the degeneracy ('density of states') of the macrostate {E,M), i.e. the number of 
configurations a with energy E and magnetisation M, and the partition function can be 
written as Z{i3, h) = J2e 9{^j ^) exp[— — /iM)]. This approach was previously 
used in the context of metastability in jHj ; see also [1^] for a method to derive suitable 
coarse-grained quantities. 

This is a useful representation, since we can compute the joint density of states 
g{E, M) numerically using the Wang-Landau algorithm fOlE!; we use a more efficient 
version of this algorithm given in JHl- (Computing g{E,M) analytically would be 
equivalent to solving the Ising model in external field, a still-unsolved problem.) The 
fact that we require the joint density of states as a function of the two parameters E 
and M restricts us to small systems [12], but we can obtain g{E, M) relatively easily for 
a system of size 32 x 32 spins, where metastability can be clearly seen under Metropolis 
dynamics. All numerical results we present are for this system size, for which the range 
of possible values for E is [-2048,2048], and for M is [-1024, 1024]. Simulation times 
are measured in Monte Carlo steps per spin (MCSS). 

From g{E, M) we can obtain the complete partition function, and hence all 
thermodynamic information at equilibrium for given values of (3 and h ^Tj. For 
example. Fig. ^ shows the shape of the equilibrium distribution —\aPQ{E,M) = 
— In g{E, M) + (3{E — hM) +lnZ for a particular f3 and h for which a metastable 
state exists. 

Two minima of different heights can be seen, separated by a saddle; the higher 
minimum corresponds to the metastable state, and the lower one to the equilibrium 
state. Fig. Q can be viewed as a 'free energy' landscape. If the system starts in the 
metastable state, then in order to escape to equilibrium, it must pass over the free 
energy barrier near the saddle point |14t 120] . 

We remark that an alternative coarse-graining has also been used to study 
metastability in the Ising model, using only the magnetisation as a coarse-grained 



Po{E,M) 



g{E, M) exp [-I3{E - hM)] 



(47) 



where 




(48) 
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Figure 1. Part of the 'free energy' ~\nPo{E,M) as a function of E and M for 
/3 — 0.5 and h = —0.02, evaluated from g{E, M) data obtained using the 2-parameter 
Wang-Landau algorithm. Two minima can be seen: the higher, metastable minimum 
is marked. Note that the z axis is logarithmic, so that there is a difference of many 
orders of magnitude between their heights. A contour plot is also shown; here the 
saddle point, the two minima, and the non-existence of certain (iJ, M) macrostates are 
visible. 

quantity. This can be motivated by considering the Ising model in the lattice gas 
representation, that is, with spin 1 representing a particle and spin —1 a void. In that 
case, the canonical ensemble is one in which M and f3 are constant, and the free energy 
is given by 

F{M; (3) = ~\nJ29iE, M)e-^^. (49) 

^ E 

Returning to the Ising model, if we now impose a magnetic field h, then the 
corresponding free energy becomes F{M; h) = F{M; (3) — hM. This can be obtained 
from the distribution of Figure ^ by summing over all E\ it is plotted in Fig. |21 F{M) 
is proportional to the logarithm of the distribution of the order parameter M fl^ |^ 
and can be calculated using several Monte Carlo methods [221 122] ■ The Wang-Landau 
method again has the advantage that we can calculate F(M; (3, h) for any parameters 
j3 and from a single run. 

We see that the free energy F{M) is still significantly non-convex. This does not. 
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Figure 2. Free energy per spin F{M;(],h)/N as a function of magnetisation per 
site M/N for L = 32, (3 = 0.5 and several values of h, obtained using the 2- 
parameter Wang-Landau algorithm. Again a metastable and an equilibrium minimum 
are visible; the former disappears for sufficiently large \h\. The inset is a comparison of 
F{M; (3 = 0.5, h = 0.05)/A^ for different system sizesL = 16, 20, 32 again as a function 
of M/N, showing the convergence towards a convex function as the thermodynamic 
limit is approached {L — > oo). 



of course, contradict the rigorous results of which show that the free energy per spin 
must be convex in the thermodynamic hmit. Indeed, the inset of Fig. |21 illustrates how 
this convexity is approached as the system size L increases. However, it shows that our 
simplified description cannot hold for arbitrarily large systems. As mentioned previously, 
to describe the metastable region adequately, we need to use macrostates specific enough 
to decide whether a critical droplet is present or not. What we are suggesting, however, 
is that a finite system described in this fashion will display significant non-convexities in 
the free energy as defined here, since there will always be a local minimum corresponding 
to the metastable state Qi{o'). 

4-2. Finding the metastable region and calculating C{a) 

We are interested in the structure of metastable states. According to our formalism, 
such states, when they exist, should be described by Qi(cr) = ZiPo^a) for configurations 
a in the metastable region T^, with C, and hence Zi, being constant over this region. To 
test this, we again look at the coarse-grained level, summing over the {E, M) macrostate 
to give 

Q,{E, M) = Z,{E, M)Po{E, M) = [1 + C{E, Mf]Po{E, M), (50) 
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where C{E,M) and Zi{E,M) are the mean values of C{a) and Zi{(t), respectively, for 
a in the {E, M) macrostate, and Qi{E, M) is the sum of Qi(cr) for a in that macrostate. 
Taking logarithms and using ln[l + C{E,M)'^] ~ 2\nC{E,M) for C{E,M) large, we 
obtain 

\nC{E,M) = ^[\nQi{E, M) - In g{E, M) + \n Z + P{E - hM)] . (51) 

If the theory is correct and, furthermore, if the parameters E and M provide an 
adequate representation of the complete phase space of the system we are studying, then 
for an {E, M) macrostate whose configurations a are all in the metastable region T^, 
we expect that C(cr) = C is constant over the macrostate, so that C{E,M) = C. We 
thus expect to have a large region in a plot of C{E, M) where it is essentially constant, 
i.e. a plateau. This region of {E, M) space, which we denote by T^, then corresponds 
to the metastable region Fm in the complete phase space. 

To find this metastable region in the reduced parameter space (for given values 
of (3 and h), we must obtain the metastable probability distribution Qi{E, M), i.e. 
the probability that the system is in the [E, M) macrostate while it remains in the 
metastable state. To do so, we record a histogram of the number of visits to each 
{E, M) pair while the system remains in the metastable state, averaging over different 
runs if necessary. Normalising this histogram then gives an estimate of the probability 
distribution Qi{E, M). It is very strongly peaked in a small region of the {E, M) plane: 
for example, for the parameters used in Fig. El the maximum value of Qi occurs at 
(Eo, Mo) = (-2040, 1022), and is given by (5i(^o, Mo) = 0.218, so that the system is in 
this single macrostate for nearly a quarter of the time spent in the metastable state; and 
adding another two macrostates gives more than half the total probability. Intuitively, 
the metastable region Fm should consist of those {E, M) pairs which have an appreciable 
Qi probability. 

We now use the Wang-Landau algorithm to calculate the joint density of states 
g{E, M) and the partition function Z for the same lattice for which we calculated 
Qi{E,M), and substitute these values into (j5T|l to obtain In C(£^,M) as a function of 
E and M. Note that this key application of the Wang-Landau algorithm determines E 
and M as the macroscopic variables to be used. 

Fig. El shows a plot of In \ C{E,M)\ for values of (3 and h such that no nucleation 
event occurred during the (long) simulation, so that the system was always in the 
metastable state. In confirmation of the theory, a large plateau is apparent. For some 
{E,M) macrostates, C{E,M) is larger than this plateau value. This happens, even 
though according to the theory it cannot since the plateau value of C is its largest 
possible value, due to the fact that these macrostates are visited very rarely during the 
simulation, so that good statistics cannot be acquired, and their measured frequency 
is larger than their true frequency. Outside the metastable region accessible in the 
simulation, we plot — InC for comparison, since there we expect that Qi{E, M) = —1/C 
(see (j24l) ). This neglects the boundary region between the two phases, to which we have 
no access using this method. In the next subsection we present an alternative approach. 
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Figure 3. \n\C{E, M)\ as a function of E and M near tlie metastable region for 
/3 — 0.8 (i.e. T ~ 0.55Tc) and = —0.1, calculated using (|5T|l . from a single run of 
7 X 10* MCSS with no nucleation events. Outside the metastable region, — In C is 
shown for comparison. 



Furthermore, the plateau value of C should be expressible in terms of the 
equilibrium distribution Pq as (c.f. (j211)) 



InC 



(52) 



^ lnPo(^,M)- Yl lnPo(^,M) 

UE,M)^I\n iE,M)ef\n 

which can be interpreted as the difference in free energy between the equilibrium and 
metastable phases. Indeed, for the case shown in Fig. 01 the plateau value calculated 
from the metastable distribution Qi using at {Eq, Mq) (where the statistics are best) 
is \nC{EQ, Mq) = 81.595, whereas the free energy difference ()52j) gives InC = 81.591. 
Note that if we write InC as a free energy difference, then it is entirely determined by 
the equilibrium distribution. The effect of the dynamics is hidden in the determination 
of the metastable region T^. 

We remark that for higher temperatures, the system does escape from the 
metastable state during a run. In this case, the identification of the metastable region Fm 
is less obvious. We take it as being those {E, M) with C{E, M) within ±1 of C{Eq, Mq). 

These results provide numerical confirmation that the metastable distribution Qi 
is proportional to the equilibrium distribution Pq in the metastable region, and that the 
proportionality constant C can be related to the difference in free energy between the 
two phases. 
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4-3. Structure ofC{a) 

To gain more insight into the function C(cr), we can use the result (j32|) . which shows that 
if we start from an initial configuration do such that C(cro) = pC, then the probability 
that after a short relaxation time the system is in the metastable state is p, while the 
probability that it is in equilibrium is 1 —p. We cannot calculate values of C(cr) directly, 
other than in the metastable region Fm, but we can use this result 'in reverse' to obtain 
a coarse-grained picture of C{a), as follows. 

For each {E, M) macrostate, we wish to generate configurations a lying in that 
macrostate. This is non-trivial, but can be accomplished by starting from a random 
initial configuration, with each spin being up or down with probability 1/2. From 
there we propose random spin flips, accepting only those which move us towards the 
desired value of {E,M). This process may get stuck, however, before reaching {E,M), 
in which case we employ Wang-Landau sampling (which is known to explore parameter 
space reasonably efficiently ^Zj) in a window containing the current and desired {E, M) 
values, to force the system into a configuration belonging to the required macrostate. 
If this does not succeed after a certain number of steps (for example if there are no 
configurations in the target 'macrostate'), then we continue with the next macrostate. 
We cannot guarantee that this procedure samples initial configurations within {E, M) 
uniformly, but empirically this seems to be the case, with no particular bias in the 
procedure. 

We start with uq configurations within the {E, M) macrostate as above, evolve each 
under Metropolis dynamics for a short time, and record whether the system has reached 
the equilibrium state, taken to be configurations with M(o") < 0, or not. The ratio 
neq/no of the number of times equilibrium is reached to the total number of trials is an 
approximation to 1 — p{E, M) for that macrostate. Note, however, that the fact that 
we average over macrostates means that we may not correctly identify the boundaries 
of the metastable region: a single macrostate may contain some configurations which 
always lead to the metastable state, and others which always lead to equilibrium, for 
example. Nonetheless it gives a clear picture of the metastable and equilibrium regions, 
and an idea of the structure of the boundary between them. 

Fig. m shows p{E, M) calculated in this way. We see clear metastable (p = 1) and 
equilibrium {p = 0) regions, separated by a boundary region where p takes intermediate 
values. The boundary region is larger than we might expect, due to the smoothing 
described above, but the system spends little time in this transition region when the 
dynamics is taken into account. Note, however, that according to the results in Sectional 
exactly where we impose the boundary between the metastable and equilibrium regions 
does not affect the results. 

Also shown in the figure is the metastable region obtained in Monte Carlo 
simulations, as described in the previous subsection. Note that the region of {E, M) with 
p{E, M) < 1 is significantly larger than this latter definition of the metastable region. 
This reflects the fact that there are configurations cr which are never reached from an 
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Figure 4. Probability p{E, M) of reaching the metastable state starting from the 
{E,M) macrostate, for (3 = 0.6, h = —0.1, and no = 50 trials for each {E,M). The 
crosses at bottom right indicate the extent of the metastable region obtained from 
Monte Carlo simulations, defined as those {E,M) having \nC{E,M) within ±a of 
lnC{Eo, Mq), with a = 1. Changing the tolerance a in the definition changes the 
horizontal extent of this region. 



initial configuration with all spins up, since the probability of doing so is negligible, and 
yet which will decay into the metastable state if started there, thus belonging to the 
metastable region according to our definition. It should also be noted that the boundary 
of the region from the simulations lies at approximately p{E, M) = 0.5, and does not 
significantly vary if the exact definition of the region is changed, in accordance with the 
results of Section El 

4-. 4- Relation of C to hysteresis loops 

Since In C corresponds to a difference in free energies, differentiating it with respect to 
the external field h gives a difference in magnetisations between the two regions: 

(53) 



dh 2 V'-'^^ 
where (M)^^ denotes the mean magnetisation in equilibrium, given by 

(^)e. := ^ (54) 

and (M)^ is similarly the mean magnetisation in the metastable state. The quantity 
(M)j^ — (M)gq has a physical meaning for those values of the external field h for which 
a metastable state exists, namely the distance on an averaged hysteresis loop between 
the two branches. 

We evaluate C as a function of h in three different ways and take the numerical 
derivative. The first is C{Eq,Mq) (i.e. C{E,M) evaluated where the metastable 
distribution Qi attains its maximum). The other evaluations use the free energy 
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Figure 5. Comparison of hysteresis loop and derivative of In C with respect to h, for 
f3 = 0.55. For each value of h, 10^ MCSS were used to find the metastable distribution. 
The hysteresis loop was obtained by averaging 1000 runs, in each run increasing h in 
steps of 0.002 and allowing the system to equilibrate at the new value of h for 5 MCSS. 
Shown are the data for the three different ways of calculating C referred to in the text: 
despite the use of a numerical derivative, the maximum deviation of these from the 
hysteresis loop data is of the order of only 2%. The inset shows the complete hysteresis 
loop and the difference between the two branches, as well as the data of type (i). 



difference dSSJ, taking to be (i) those {E,M) for wliich C{E,M) is within ±1 
of C{Eq, Mo), and (ii) those {E,M) for which M lies on the right of the maximum 
of F{M;/3,h), which is a more 'traditional' method |2a|. Fig. El shows {2/(3){dC /dh) 
compared to the difference between the heights of the two branches on a hysteresis loop. 
The agreement is reasonably good, including for the more traditional method, although 
the data from C{Eq, Mq) is noisy. We remark that this provides a possible experimental 
avenue for measuring a physical quantity directly related to C. 

5. Several metastable states 

In contrast to the systems we have discussed thus far, it may happen that a given system 
has several metastable states: see for example studies of a Blume-Capel model with 
two metastable states in and also (2^1 • In this section we extend our formalism 

to describe such situations. As before, instead of focusing on a specific example, we 
approach the problem through the general formalism of Markov processes satisfying 
detailed balance. Previous results in a similar direction can be found in Refs. |^ I27j. 
In the following, we limit ourselves to the case in which the number of metastable 



Metastability in Markov processes 



23 



states is independent of the system size A^. Other situations are also possible: for 
example, it is generally assumed that the physics of both structural and spin glasses may 
be related to the presence of a macroscopic number of metastable states j2H!- However, 
such a scenario presents significant additional complexities which we do not address. 
In particular, it is not clear that for such systems there really exists an appropriate 
description in thermodynamic terms, as we show in this paper for the systems we call 
metastable. 

Since the following may well appear unnecessarily complex, let us first explain the 
origin of the difficulties that may arise when dealing with multiple metastable states. In 
the case in which only one metastable state is present, there is only one eigenstate Pi, 
which is essentially non-zero in the metastable region. To generalize this to the case of K 
isolated metastable states, all of which decay to equilibrium, is indeed straightforward: 
one then finds K different regions and K eigenstates, one concentrated on each region, 
and everything is essentially very similar to the case of a single metastable state. The 
non-trivial issue arises when one metastable state must nucleate another metastable 
state before it can reach equilibrium. Under these circumstances, there is no clear 
correspondence between the regions in which the Pa are significantly different from zero 
and the metastable regions. We must therefore proceed slightly differently, as follows. 

We now denote by all the eigenstates of the operator L of the master equation 
(jni which have small relaxation rates Qa- The various Qa may either be all of the same 
order, or differ considerably from each other. The crucial point is that they satisfy 
Qq. -C i-e. they are all "small", and their number K should be fixed, independent 

of system size A^. 

In analogy to the case of systems with a single metastable state, we define 
Pa{cr) 



CJa] 



\Ca\ '■= max |Ca(cr)| , 
o-Gr 

r„ := {cj:{1- K)Ca < Ca{cj) < C^} . 

Note that the eigenstates Pa{cr) are defined up to a global sign; we choose the sign so 
that the Ca are positive. The numbers < < 1 are chosen to ensure that the sets 
Fq are disjoint. For these states to be metastable, we will assume that such a set of 
numbers exists, and that they are 0(1). 

Using exactly the same approach as in the previous section we can show that 
e^"*CQ,[cr(t)] is a martingale for any a. Defining as the first time that the system 
leaves F^, given that it starts with an initial condition (Tq, such that Ca{<ya) = Ca, then, 
as before, 

P(T„ < t) < ^ (1 -e~^"*) = O(fiat). (55) 

Aq 

Thus, if we consider the initial distribution ^^-^(o"), then, after an equilibration time 
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of order the system will be described by a probability distribution given by 

K 

Qa{(r) ■.= Po{a) + J2Cp{a^)Pp{a). (56) 

/3=1 

Due to (jSSj), the probability that the process beginning at cTq, leaves Fq, in the 
relevant time range is very small, so that 

Being a probability distribution, Qa{o') is non- negative, so the above result implies that 
Qo(c") ~ for cr ^ Fq,. Thus, since the regions Fq are disjoint, we conclude that for all 
a e F, 

Qa{cr)Qp{a) ^0 for « ^ /3 (58) 

and also that each is normalized over the region Fq,. These functions play the role 
of Qi in the case with a single metastable state. 

Again, it is straightforward to show that a restricted process can be constructed 
inside each Fq,, and that such a process remains close to the original unrestricted process 
in the relevant time range if the initial condition of both processes satisfies Ca{cra) = Ca- 
Thus, we can identify 

= ZaPoi^r) (a G F,), (59) 



where the constant Za is given by 



(60) 



in analogy to equation (j24|) . We can therefore again interpret as the partition function 
of the ensemble restricted to F^. 

Defined in this way, the are orthogonal to each other although, being normalized 
as probability distributions, they are not orthogonal to Pq. Nevertheless, the functions 
Qa together with Pq still form a linearly independent basis set, and the description of the 
system can be carried out in terms of these functions, which are essentially stationary 
in the relevant time range. 

Now, for the slowly decaying states Qa to describe metastability, an additional 
condition is still required, namely 

« 1 (61) 



for all a. Thus, the functions Qa are not only assumed to essentially be different from 
zero on disjoint sets, but also, the states they describe are assumed to be extremely 
improbable in equilibrium. 

The aim now is to show that, if the system starts from an arbitrary initial condition, 
then, with high probability, it either evolves to a state for which the Caicr) are close to 
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those of a Qa{<^) or to equilibrium. Further, this happens on a "short" timescale, that 
is, of the order of at most 

As before, the first step in this direction is to consider the evolution of a system 
starting at an arbitrary a^. After a time of order has elapsed, the system will be 

described by the probability distribution 

K 

P{a\ao) = Po{ct) + Ca(ao)P„(a). (62) 

a=l 

We can express the functions in terms of the Qa(cr) and the equilibrium 

distribution Pq{<7), as 

= E - Poia) f: (63) 

/3=1 13=1 

where the coefficients are obtained from (jSUj) and ()59|) and from the orthogonality 
between Po{<y) and Paicr), as well as from the fact that -Po(c") is negligible on each 
Fq. Thus we can rewrite the expression for P{a\ao) as 



P{a\ao) 



f3=l a=l 



Caicro)Caio-f3) 



Qp{a). (64) 



13=1 \_a=l 

The above expression means that the system has evolved to one of the states described 
by a Q/3 distribution, with probability 

K 



P{a e F^lao) = E ^(^1 



a=l 



(65) 



and will decay to equilibrium with probability 



P(^GreJcro) = l-EE 

(3=1 a=l 



(66) 



Thus, for the process a{t) starting at cr(t = 0) = ctq, the functions Ca[(T(t)] tend to the 
values Caicp) with probability P(cr G T/3\ao), or to zero with probability P{a G Feq|(Jo), 



in a time of order 



Further, note that if we choose ao = cr^ in equation (jHS)), then we find 



Ca{cr-y)Ca{crf3) 



a=l 1^ 



6. 



7/3- 



(67) 



Otherwise, the system would quickly leave F^ even though it had started at a^, which is 
contrary to what we have shown in equations ()55|) and (jUT|) . This can also be understood 

— 1/2 

as a statement that the functions Za Qa are othonormal with respect to the scalar 
product (jD). 

Thus, a system starting at an arbitrary initial state will either decay to equilibrium 
or evolve to a state described by one of the metastable distributions Qa- If it reaches the 
metastable state Qa, then the function Ca{cr) for this process will grow to Ca = Ca{cra)- 
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This occurs in a time of order after which the results pertaining to processes 

characterized by having the maximal value of Ca{<y) apply- In particular, it will be very 
probable that the process remains in Fq, for a long time. Further, once in the metastable 
state, it is described by the equilibrium distribution restricted to that region of phase 
space, which is again the expected physical behavior of metastable states. 

Finally, note that a formula analogous to Zi = 1 + C^, which was derived for the 
case of one metastable state, can be obtained as follows. We have 



(68) 



1 + J]C7^(ajC^(a) 

and upon substituting cTq, for a and using the fact that Qa{<7a) = ZaPo{<Ja), one finally 
obtains 

Z^ = l + J2Cf3iaaf. (69) 

The other feature that can be discussed within this formalism is the decay path of 
a metastable system. For this concept to be clear cut, we will consider the case in which 
the decay rates of the metastable states, while still small, are very different from each 
other. In such a scenario, there are time scales on which the fastest metastable state 
has decayed with certainty whereas no other one has. The question then is whether we 
can evaluate the probability that a system originally in the short-lived metastable state 
will decay to another metastable state. 

Since the eigenvalues of the evolution operator are assumed to be ordered in 
increasing magnitude, we denote by Qk{o') the distribution describing the metastable 
state with shortest lifetime. This distribution has support on F^^ which is characterized 

Now, after a time greater than ^j^^, the contributions from Pk{<^) vanish and the 
initial distribution evolves into 

Q'M = Po(a) + Yl CpMPpia), (70) 

f3=l 

with the probability of finding the system in the initial metastable state vanishing: 

although Q'j^ is still normalized. In particular, this means that Qk{o-) = C'x(o"i<')-PA'(<7) 
for a G Tk- 

Expressing Q'xi^') in terms of the remaining Q^((t), we find that the probability 
that the state be found in Fq, is given by 

p Ci3{crK)Cp{cro,) CK{crK)CKM 

r^K^a - > ^ ^ , { I ^) 

where we have used (jFTfjl in the last equality. 
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Thus, the probabihty for the decay from this metastable state to another can be 
obtained from the values of the coefficients Caicr^) appearing in the previous equation. 
(Note that the above expression imphes that CxicTa) < 0.) The probabihty that the 
state instead decays directly to equilibrium is 



Thus, under certain circumstances, when the decay rates of the various eigenmodes are 
very different or their corresponding metastable states fall along different paths, the 
complete decay of the system can be described as an irreversible Markov chain with 
transition probabilities given by equations (f7^ and (fTSj) . 

6. Conclusions and Outlook 

A fundamental issue in metastable systems concerns the possibility of describing them 
as thermodynamic equilibria (in an extended sense) of the system at hand. The fact that 
they decay to an equilibrium state which is in general quite different appears to preclude 
such an approach, as do the various results concerning the impossibility of analytically 
continuing the free energy beyond the coexistence curve. On the other hand, the use of 
a thermodynamic approach is routine in the applied work on the subject. 

In this work, we have attempted to justify the thermodynamic approach starting 
from a Markovian description. While this may, at first, seem to be an exceedingly 
restrictive assumption, a moment's thought will show the contrary: indeed, the 
Markovian approximation is expected to become reasonable on large time scales. 
But metastability is essentially concerned with that time range which covers times 
much larger than any microscopic relaxation time but much shorter than the average 
nucleation time (which we have, throughout the paper, called "the relevant time range"). 
If such a range does not exist, or if it is not large enough, then we may not meaningfully 
speak of a metastable state. The Markovian approximation is therefore expected to be 
relevant within the time range of interest. 

Further, we assume that the Markovian dynamics satisfies detailed balance. This 
last condition is not essential, and indeed, other approaches along similar lines have 
been made without the assumption of detailed balance EII- However, in addition 
to simplifying the derivation, this assumption allows the unambiguous identification of 
the probability distribution describing the equilibrium state with that of equilibrium 
statistical mechanics, thus justifying the use of concepts from equilibrium statistical 
mechanics and thermodynamics in the description of metastable states. 

Finally, we have limited ourselves systematically to finite systems. On the one 
hand, this comes from the fact that severe technicalities arise whenever infinite systems 
are considered as such. More specifically, however, as already argued by Gaveau 
and Schulman the thermodynamic limit presents some unique difficulties for 
metastability: indeed, since the size of the critical nucleating droplet remains constant 



K-l 




(73) 



Metastability in Markov processes 



28 



as the thermodynamic hmit is attained, the time in which the first such droplet will 
arise goes to zero as the system size increases to infinity. 

Our definition of a metastable state includes two components: first, it should involve 
an isolated eigenstate (when the system has only one metastable state) of the master 
operator having an exceptionally low eigenvalue. This eigenstate allows us to define 
the metastable region in the phase space of the system. Second, we impose a technical 
condition (fT^ meaning that the probability of being in a metastable state at equilibrium 
is vanishingly small. The first condition serves to discard the possibility of slow decay 
mediated by long wavelength hydrodynamic modes. Indeed these usually arise as a 
quasi- continuum of low-lying excitations and therefore cannot satisfy the condition of 
being isolated. On the other hand, the technical condition (fT^ is rather more difficult to 
understand. Physically, however, since the properties of the metastable and equilibrium 
phases are markedly different, it is clear that we should make such a requirement of any 
metastable state, as was already pointed out in jl]. 

Under the above hypotheses we have shown the following results: first, that any 
initial condition whatsoever will relax, in a short time, to a state which is either fully 
in the metastable region or to equilibrium. Further, any state starting well inside the 
metastable region has a very low probability of leaving it in the relevant time range. 
We were also able to generalize these results to the case in which a finite number of 
metastable states exist. We could not, however, extend this to situations in which the 
number of metastable states grows with the system size: this clearly cannot be done, 
since it would include, among others, the case of slowly decaying hydrodynamic modes, 
which correspond to a physically entirely different situation. 

A further important result allows to justify the thermodynamic treatment: we show 
that if one starts inside the metastable region, then a Markov process which reflects the 
system whenever it attempts to leave the metastable region, is in fact close to the 
original physical process over the relevant time range. We may therefore, for properties 
which can be observed over the relevant time range, use this "restricted" process instead 
of the original one: all the difficulties associated with the existence of nucleation then 
disappear, so we may apply the entire machinery of equilibrium statistical mechanics 
to it (Green-Kubo formulae, linear response and so on) while remaining close to the 
correct answer for the original system. 

The main open issue clearly concerns systems with a macroscopic number of low- 
lying eigenstates of the master operator. At least two apparently different classes of 
such systems are known: on the one hand, as we have said before, systems in which 
slow hydrodynamic modes play a role. On the other hand, both structural and spin 
glasses are assumed to exhibit a large number of metastable states. Clearly, neither 
can, at present, be treated by the methods presented here, but their extension to such 
systems certainly presents an interesting challenge. 
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